设计矩阵 in R

1.利用设计矩阵创建复杂模型

问题引入1：如前面所述，我们想了解在考虑小鼠体重的情况下，突变组与对照组的小鼠体积是否有差异。我们使用设计矩阵与方程结合得到方法对数据进行拟合，并探讨方程的F值和p值，可参考往期推文设计矩阵（design matrices）。

#### 利用设计矩阵，实现t检验与回归联合分析 #######
# 1.创建对照组和突变组的标签（labels）
Type <- factor(c(
  rep("Control", times=4),
  rep("Mutant", times=4)))
# 2. 依次添加对照组和突变组的数据（weight与size相对应，代表来自同一只小鼠的数据）
Weight <- c(2.4, 3.5, 4.4, 4.9, 1.7, 2.8, 3.2, 3.9)
Size <- c(1.9, 3, 2.9, 3.7, 2.8, 3.3, 3.9, 4.8)

# 3.查看设计矩阵（标准设计矩阵）
model.matrix(~Type+Weight)  

#注释：
#波浪号代表y轴的数据（这里指的是小鼠的体积）；
#小鼠的体积由对照组均值（control(mean)），小鼠的类型（Typemutant）和小鼠的体重(weighgt)创建模型 ；
#在设计矩阵中，R默认对照组均值control(mean)参数（对应于设计矩阵的第一列）。

#如下为设计矩阵的详细信息：
#   (Intercept) TypeMutant Weight
# 1           1          0    2.4
# 2           1          0    3.5
# 3           1          0    4.4
# 4           1          0    4.9
# 5           1          1    1.7
# 6           1          1    2.8
# 7           1          1    3.2
# 8           1          1    3.9
# attr(,"assign")
# [1] 0 1 2
# attr(,"contrasts")
# attr(,"contrasts")$Type
# [1] "contr.treatment"

#设计矩阵中的第一列乘以control(mean)，第一列均为1，代表该列在所有数据中均为开（on）；
#设计矩阵中的第二列乘以突变组与对照组之差，1代表数据在突变组中为开（on），表示突变组的数据；0对应对照组数据，表示突变组的数据；
#设计矩阵中的第三列乘以slope，代表考虑小鼠体重与小鼠体积的关系。

# 4.lm（）进行回归分析
model <- lm(Size~Type + Weight)

# 查看模型的统计量
summary(model)  

# Call:
#   lm(formula = Size ~ Type + Weight)
# 
# Residuals:
#   1        2        3        4        5        6        7        8 
# 0.05455  0.34562 -0.41623  0.01607 -0.01753 -0.32646 -0.02062  0.36461 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)   
# (Intercept)  0.08052    0.52744   0.153  0.88463     #对单个变量的检验
# TypeMutant   1.48685    0.26023   5.714  0.00230 **  #对单个变量的检验
# Weight       0.73539    0.13194   5.574  0.00256 **  #对单个变量的检验
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.3275 on 5 degrees of freedom
# Multiple R-squared:  0.8975, Adjusted R-squared:  0.8564 
# F-statistic: 21.88 on 2 and 5 DF,  p-value: 0.003367  #对模型的整体检验

模型统计量的解读

• 对模型整体进行F检验， 其相对于仅考虑小鼠体积的均值拟合直线，F值为21.88，p值= 0.003367，说明考虑小鼠体积与小鼠体重关系、考虑小鼠类型的模型的预测效果好于仅考虑小鼠体积均值拟合直线的预测效果。

• 对模型中weight参数的解读： p=0.00256，说明考虑小鼠体重与小鼠体积关系的模型，较不考虑小鼠体重与小鼠体积关系的模型的预测效果好。

• 对模型中TypeMutant参数的解读： p=0.00230，说明考虑小鼠类型的模型，较不考虑小鼠类型的模型的预测效果好。

结合回归分析的统计检验结果，我们可以对最初的问题作出回答：同时考虑小鼠体积与小鼠体重关系、小鼠类型的模型的预测效果更好，故考虑小鼠体重与小鼠体积的关系时，突变组与对照组的小鼠体积具有显著差异。

2.利用设计矩阵实现批次校正

问题引入2：在考虑批次效应时，整个模型的检验结果并不是我们关注的结果。例：由于批次效应，lab B在重复lab A的工作时，发现其总体测量结果较lab的总体测量结果偏低。如果要结合这两个数据集，分析是否突变组数据与对照组数据有显著差异，我们首先需要校正批次效应。

# 1.创建lab和type标签
Lab <- factor(c(
  rep("A", times=6),
  rep("B", times=6)))
Type <- factor(c(
  rep("Control", times=3),
  rep("Mutant", times=3),
  rep("Control", times=3),
  rep("Mutant", times=3)))

# 2.依次添加数据
ExpreSSion <-c(
  1.7, 2, 2.2,
  3.1, 3.6, 3.9,
  0.9, 1.2, 1.9,
  1.8, 2.2, 2.9)

# 3.查看设计矩阵（用于明确设计矩阵是否正确）
model.matrix(~Lab+Type)
#    (Intercept) LabB TypeMutant
# 1            1    0          0
# 2            1    0          0
# 3            1    0          0
# 4            1    0          1
# 5            1    0          1
# 6            1    0          1
# 7            1    1          0
# 8            1    1          0
# 9            1    1          0
# 10           1    1          1
# 11           1    1          1
# 12           1    1          1
# attr(,"assign")
# [1] 0 1 2
# attr(,"contrasts")
# attr(,"contrasts")$Lab
# [1] "contr.treatment"
# 
# attr(,"contrasts")$Type
# [1] "contr.treatment"

# 4.lm()回归分析
batch.lm <- lm(ExpreSSion ~ Lab + Type)

# 查看回归分析的结果
summary(batch.lm)
# Call:
#   lm(formula = ExpreSSion ~ Lab + Type)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -0.6500 -0.2833 -0.0500  0.2750  0.7167 
# 
# Coefficients:
#           Estimate Std. Error t value Pr(>|t|)    
#   (Intercept)   2.1167     0.2279   9.287  6.6e-06 *** #对单个变量的检验
#   LabB         -0.9333     0.2632  -3.546 0.006250 **  #对单个变量的检验
#   TypeMutant    1.2667     0.2632   4.813 0.000956 *** #对单个变量的检验
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4558 on 9 degrees of freedom
# Multiple R-squared:  0.7989, Adjusted R-squared:  0.7542 
# F-statistic: 17.87 on 2 and 9 DF,  p-value: 0.0007342  #对整个模型的检验

模型统计量的解读：

• 对模型的整体F检验： 在这个问题中，我们关心的不是整个模型的预测效果，而是校正批次效应后，突变组和对照组的基因表达有无显著差异。故最后F检验的结果不是我们的关注对象。在此结果中，p值=0.0007342，说明两个实验室之间是存在批次效应的。

• 对TypeMutant参数的检验： p=0.000956，提示含有突变类型（TypeMutant）的模型的预测效果优于不含突变类型模型的预测结果，即考虑批次效应和突变类型的模型的预测效果优于仅考虑批次效应的预测效果。基于此，我们可以得出结论：在校正批次效应后，突变组和对照组的基因表达量存在显著差异。

参考视频：

https://www.youtube.com/watch?v=Hrr2anyK_5s&list=PLblh5JKOoLUIzaEkCLIUxQFjPIlapw8nU&index=8

编辑：吕琼

校审：罗鹏